#ifndef DRAW_H
#define DRAW_H

const double PI = 3.1415926 ;
const unsigned WINDOW_WIDTH = 500;
const unsigned WINDOW_HEIGHT = 600 ; 

const float gaussian_sigma = 0.84089642f ;
const float gaussian_constant[5] = 
//{0.22508352f,	0.11098164f, 0.01330373f, 0.00038771f, 0.0f};
{0.22508352f,	0.11098164f, 0.05330373f, 0.0138771f, 0.0f};



float randreal() ;


//void prepare_color(std::vector<Face *> & faces); 
struct int_RGB
{
	int_RGB (unsigned _r, unsigned _g, unsigned _b): r(_r), g(_g), b(_b) {}
	int_RGB():r(255),g(255),b(255) {}
	unsigned r, g, b;
	friend inline int_RGB operator * (const int_RGB & color, float scale)
	{
		return int_RGB(int(color.r * scale), int(color.g * scale), int(color.b * scale)) ;
	}
	friend inline int_RGB operator + (const int_RGB & lhs, const int_RGB & rhs)
	{
		return int_RGB(lhs.r + rhs.r, lhs.g + rhs.g, lhs.b + rhs.b);
	}
};

inline float linear_gaussian_blur(float distance_ratio) /*[0,1] from the center*/
{
	int a = int(distance_ratio * 4);
	int b = a + 1;
	if (b >= 5) b = a;
	float weight_a = distance_ratio * 4 - a;
	float weight_b = 1.0f - weight_a;
	return ( gaussian_constant[a] * weight_b + gaussian_constant[b] * weight_a ) / 0.22508352f  ;
}
struct color_RGB
{
	color_RGB(float _r, float _g, float _b):r(_r),g(_g),b(_b) {}
	color_RGB():r(1.0),g(1.0),b(1.0) {}
	float r, g, b;
};
int_RGB rand_color()  ;

float compute_texture() ;
struct STATES
{
	STATES()
	{
		sample_points = false;
		selected_face_idx = 0; 
	}
	bool sample_points ; /*Controled by key 'S'*/
	unsigned selected_face_idx ;
};
class Point2 {
public:
	Point2(float x_in, float y_in) : x(x_in), y(y_in) { }
	Point2() : x(0), y(0) { }
	float x ;
	float y ;

	inline float length() const
	{
		return sqrt (x*x + y * y);
	}

	friend inline std::ostream& operator<<(std::ostream& out, const Point2& v) 
	{
		return out << v.x << " " << v.y ;
	}

	friend inline std::istream& operator>>(std::istream& in, Point2& v) 
	{
		return in >> v.x >> v.y ;
	}

	friend inline Point2 operator+(const Point2 & v1, const Point2 & v2) 
	{
		return Point2( v1.x + v2.x, v1.y + v2.y) ;
	}

	friend inline Point2 operator-(const Point2 & v1, const Point2 & v2) 
	{
		return Point2(v1.x - v2.x, v1.y - v2.y) ;
	}

	friend inline Point2 operator * (const Point2 & lhs, float rhs)
	{
		return Point2 (lhs.x * rhs, lhs.y * rhs);
	}
	friend inline float operator ^ (const Point2 & lhs, const Point2 & rhs) 
	{
		return (lhs.x * rhs.y - lhs.y * rhs.x );
	}
} ;
class Point3 
{
public:
	Point3 () {}
	Point3 (float a, float b, float c)
	{
		x = a; y = b; z = c;
	}
	Point3 (Point2 & p2)
	{
		x = p2.x ;
		y = p2.y ; 
		z = 0.0f ;
	}
	float x, y, z;
	inline float length() const
	{
		return sqrt (x*x + y * y + z*z);
	}
	inline void normalize() 
	{
		float len = length() ;
		x /= len; 
		y /= len;
		z /= len;
	}
	friend inline float operator * (const Point3 & lhs, const Point3 & rhs)
	{
		return (lhs.x * rhs.x + lhs.y * rhs.y + lhs.z * rhs.z);
	}
	friend inline Point3 operator + (const Point3 & lhs, const Point3 & rhs)
	{
		return Point3 (lhs.x + rhs.x , lhs.y + rhs.y , lhs.z + rhs.z);
	}
	friend inline Point3 operator / (const Point3 & lhs, float rhs)
	{
		return Point3 (lhs.x / rhs, lhs.y / rhs, lhs.z / rhs);
	}
	friend inline Point3 operator * (const Point3 & lhs, float rhs)
	{
		return Point3 (lhs.x * rhs, lhs.y * rhs, lhs.z * rhs);
	}
	friend inline Point3 operator - (const Point3 & lhs, const Point3 & rhs)
	{
		return Point3 (lhs.x - rhs.x , lhs.y - rhs.y , lhs.z - rhs.z);
	}
	friend inline Point3 operator ^ (const Point3 & lhs, const Point3 & rhs) 
	{
		return Point3 (lhs.y * rhs.z - lhs.z * rhs.y, 
			(lhs.z * rhs.x - lhs.x * rhs.z),
			lhs.x * rhs.y - lhs.y * rhs.x );
	}
	friend inline std::ostream& operator<<(std::ostream& out, const Point3 & v) 
	{
		return out << v.x << " " << v.y << " " << v.z ;
	}
	friend inline std::istream& operator>>(std::istream& in, Point3 & v) 
	{
		return in >> v.x >> v.y >> v.z ;
	}

	static Point3 normal (const Point3 & p)
	{
		float length = sqrt (p.x * p.x + p.y * p.y + p.z * p.z);
		return Point3 (p.x /length, p.y /length, p.z /length);
	}

};

typedef Point3 Point ;
typedef Point3 Vector ;
typedef Point3 Vector3 ;



typedef Point2 Vector2 ;


class Resolution 
{
public:
	Resolution():reso(1) {}


private:
	unsigned reso ; 
	//std::vector<std::vector<color_RGB> > color_array; 
};

inline int barycentric2index (unsigned res,  unsigned i, unsigned j)
{
	return  (res + 1 + res + 2 - i ) * i / 2 + j; 
}
void index2barycentric (unsigned res, unsigned idx, unsigned & i, unsigned & j);


//void drawcircle(float x, float y, float z, float radius);
//void drawsphere(float x, float y, float z, float radius);

static void getCoordEnv (double modelViewMatrix[16], double projectionMatrix[16], int viewPort[4])
{
	glGetIntegerv(GL_VIEWPORT, viewPort);	
	glGetDoublev (GL_MODELVIEW_MATRIX, modelViewMatrix);
	glGetDoublev (GL_PROJECTION_MATRIX, projectionMatrix);
}
static 	Vector3 mappingWin2Obj (Vector3 winPos, double * modelViewMatrix, double * projectionMatrix, int * viewPort)
{
	double x, y, z;
	gluUnProject(winPos.x, WINDOW_HEIGHT - winPos.y, winPos.z, modelViewMatrix, projectionMatrix, viewPort, 
		& x, &y , &z);
	return Vector3 (float(x), float(y), float(z));
}
inline Vector2 mappingObj2Win (Vector2 objCoord, double * modelViewMatrix, double * projectionMatrix, int * viewPort)
{
	double winx, winy, winz ;
	gluProject(objCoord.x, objCoord.y, 0.0f, modelViewMatrix, projectionMatrix, viewPort, & winx, & winy, & winz );
	winy = WINDOW_HEIGHT - winy ;
	return Vector2 (float(winx), float(winy));
}
inline Vector3 mappingObj2Win (Vector3 objCoord, double * modelViewMatrix, double * projectionMatrix, int * viewPort)
{
	double winx, winy, winz ;
	gluProject(objCoord.x, objCoord.y, objCoord.z, modelViewMatrix, projectionMatrix, viewPort, & winx, & winy, & winz );
	winy = WINDOW_HEIGHT - winy ;
	return Vector3 (float(winx), float(winy), float(winz));
}
inline bool isInTriangle (Vector2 v[3], Vector2 & point)
{
	v[0] = v[0] - point; v[1] = v[1] - point; v[2] = v[2] - point;
	float d[3];
	d[0] = v[0]^v[1];d[1] = v[1]^v[2];d[2] = v[2]^v[0];
	return (d[0] < 0 && d[1] < 0 && d[2] < 0) ;
}


#define Matrix3  float **

inline float det3 ( float m[3][3])
{
	float det = m[0][0] *(m[2][2]*m[1][1] - m[2][1]* m[1][2]) 
			- m[1][0] * (m[2][2]*m[0][1] - m[2][1]* m[0][2])
			+ m[2][0] * (m[1][2]*m[0][1] - m[1][1]* m[0][2]);
	return det ;
}
inline void inverse3 (float m[3][3],  float invm[3][3])
{
	float rdet = 1.0f / det3(m);
	invm[0][0] =   rdet * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) ;
	invm[0][1] = - rdet * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) ;
	invm[0][2] =   rdet * (m[0][1] * m[1][2] - m[0][2] * m[1][1]) ;

	invm[1][0] = - rdet * (m[1][0] * m[2][2] - m[1][2] * m[2][0]) ;
	invm[1][1] =   rdet * (m[0][0] * m[2][2] - m[0][2] * m[2][0]) ;
	invm[1][2] = - rdet * (m[0][0] * m[1][2] - m[0][2] * m[1][0]) ;

	invm[2][0] =   rdet * (m[1][0] * m[2][1] - m[1][1] * m[2][0]) ;
	invm[2][1] = - rdet * (m[0][0] * m[2][1] - m[0][1] * m[2][0]) ;
	invm[2][2] =   rdet * (m[0][0] * m[1][1] - m[0][1] * m[1][0]) ;
}
inline Vector3 multiply_vector(float m[3][3], Vector3 & v)
{
	Vector3 result;
	result.x = m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z;
	result.y = m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z;
	result.z = m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z;
	return result;
}
#endif
